Mean-field behavior of the sandpile model below the upper critical dimension 
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We present results of large scale numerical simulations of the Bak, Tang and Wiesenfeld sandpile 
model. We analyze the critical behavior of the model in Euclidean dimensions 2 < d < 6. We 
consider a dissipative generalization of the model and study the avalanche size and duration distri- 
butions for different values of the lattice size and dissipation. We find that the scaling exponents 
in d = 4 significantly differ from mean- field predictions, thus suggesting an upper critical dimen- 
sion d c > 5. Using the relations among the dissipation rate e and the finite lattice size L, we find 
that a subset of the exponents displays mean-field values below the upper critical dimensions. This 
behavior is explained in terms of conservation laws. 

PACS numbers: 64.60.Lx, 05.40.+j, 05.70.Ln 
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Since the introduction of the concept of self-organized 
criticality (SOC) ten years ago an enormous effort 

has been devoted to the understanding of this irreversible 
dynamical phenomenon. SOC models oppose the stan- 
dard picture of critical phenomena, since their dynamics 
should generate a self-organization of the system into a 
critical state, without need for the fine tuning of external 
parameters. The paradigmatic SOC model is the sand- 
pile automaton, in which a slow external driving of sand 
particles leads to a stationary state with avalanches dis- 
tributed on all length scales |3J. Despite the apparently 
simple rules, the model shows a complicated behavior 
which is not amenable to a complete solution. 

In SOC models, the concept of "spontaneous" critical- 
ity is quite ambiguous because it has been recognized that 
criticality appears only if the driving rate is fine tuned 
to zero H-0]. The slow driving assumption implies non- 
locality in the dynamical rules of the model (pf, which 
makes a general theory of SOC problematic ||. Sev- 
eral important theoretical questions are still not resolved, 
such as the precise definition of universality classes, the 
value of the upper critical dimension, and the validity 
of fluctuation-dissipation theorems. These problems are 
also reflected in the relatively few exact results available 
in the literature 0,^]. Furthermore, these issues are also 
unclear from the numerical point of view, and only in the 
last years, earlier computational efforts ||,[l0) have been 
followed by more accurate numerical studies []Tl|-p~3|] . 

Recently, a general dynamical mean-field (MF) analy- 
sis [[| of sandpile models pointed out the similarities be- 
tween SOC models and phase transitions in systems with 
absorbing states p^ |. Criticality is analyzed in terms of 
the response function singularities and the MF critical 
exponents are calculated. This method relates bulk and 
boundary dissipation and introduces a scaling relation re- 
lating dissipation and finite-size effects. Moreover, due to 
the conservative nature of sandpiles at the critical point, 



a subset of critical exponent was predicted to display MF 
values also in low dimensions Q . This result provides an 
important test to verify the validity of the MF theory, 
and can be used as a consistency check in the numeri- 
cal analysis of several exponents characterizing sandpile 
models. 

Here, we study the critical behavior of the avalanche 
size and duration distribution in order to provide numer- 
ical evidences for the MF behavior of low dimensional 
sandpiles. We perform an accurate study of critical ex- 
ponents for conservative and dissipative |l5|,[T(| sand- 
piles in dimensionality ranging from d = 2 to d = 6. 
This allows us to estimate the upper critical dimension 
d c . In contrast with recent numerical simulations [ p^[ , 
MF behavior is observed only in d = 6 and we therefore 
exclude that d c = 4. In addition we found that some 
critical exponents assume constantly their MF values in 
all dimensions d, as predicted in Ref. 0]. 

We consider the d-dimensional Bak, Tang and Wiesen- 
feld (BTW) sandpile model on a hypercubic lattice of 
size L. On each site i of the lattice we define an integer 
variable Zi which is identified with the sand or energy 
stored in the site. At each time step an energy grain is 
added on a randomly chosen site (zj — » 2, + 1). When 
one of the sites reaches or exceed the threshold z c = 2d a 
dynamical process occurs: Zi = zi — 2d and Zj = Zj + 1, 
where j represents the nearest neighbor sites. Such a 
"toppling" event can induce nearest neighbor sites to top- 
ple on their turn and so on, until all sites are below the 
critical threshold. This process is called an avalanche. 
The slow driving condition is implemented by stopping 
the random energy addition during the avalanche spread- 
ing. This means that the driving time scale is infinitely 
slow with respect to the avalanche characteristic time. 

The model is locally conservative; no energy grains are 
lost during the toppling event. The only dissipation oc- 
curs at the boundary, from which energy can leave the 
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system. We also use a nonconservative definition of the 
model. With probability p the toppling site loose its 
energy without transferring it to its nearest neighbors. 
This means that on average a quantity e = 2dp of en- 
ergy is dissipated in each toppling. In this case periodic 
boundary conditions can be considered. With both these 
definitions, the model reaches a stationary state in which 
the energy introduced by the external random drive is 
balanced on average by the energy dissipated in the dy- 
namical evolution. In the stationary state, we can de- 
fine the probability that the addition of a single grain 
is followed by an avalanche of s relaxation events. In 
the limit e — > 0, it is possible to show that the system 
response function is diverging, revealing the presence a 
critical point ||. Close to criticality, the avalanche size 
distribution assumes the scaling form 

P(s) = S - T G(s/s c ) , (1) 

where s c is the cutoff in the avalanche size. 

In the infinite time scale separation, the cut-off size 
is a function s c ~ of the bulk or border dissipa- 

tion. The boundary dissipation follows the scaling form 
e ~ where \x is the exponent relating the dissipation 
rate with the system size. Thus we obtain that in the case 
of a fully conservative system s c ~ L^/ a . It is useful to 
introduce also the avalanche characteristic length £ and 
the scaling relations s c ~ £ and £ ~ e _iy , which define 
the the fractal dimension and the characteristic length 
divergence exponents, respectively. By noting that £ and 
L must rescale in the same way, we immediately obtain 
the scaling relations: 

Da = v~ l , u = fi~ 1 . (2) 

The MF theory gives tmf = 3/2, gmf = 1/2 and 
Dmf = 4 f|]. In addition, the theory of Ref. j| pre- 
dicts that /x = 2 and v — 1/2 in all dimensions because 
of the inherent conservation law of these models. The 
values of these two exponents also imply that < s >^ L 
and < s >^ e~ 7 with 7 = 1 for any d fT^ ]. From these 
results, we obtain the scaling relation Da = 2, that also 
holds for all d. These results provide a powerful consis- 
tency check in the numerical analysis of several exponents 
characterizing sandpile models. The value of the expo- 
nents r, a and D depend on d and will only agree with 
MF theory values when d > d c . 

In order to test the above picture we have studied 
the avalanche size distribution in systems with dimen- 
sion ranging from d = 2 to d = 6 and varying sizes L and 
dissipation e. In the first simulation set (e = 0), system 
sizes L < 1024 for d = 2, L < 762 for d = 3, L < 144 
for d = 4, L < 53 for d = 5 and L < 27 for d — 6 have 
been investigated. In the second set the dissipation rates 
change with the dimension: e > 10 -5 for d = 2, e > 10~ 4 
for d — 3 and e > 10 _1 for d = 4, 5 and 6 with lattice 
of the maximum size available. In each case, statistical 
distributions are obtained averaging over a number rang- 
ing from 10 6 to 10 7 nonzero avalanches. For d > 3, the 



sizes reached in our simulations are, to our knowledge, 
the largest which have ever been used. In d = 2 we did 
not push the computational effort too far, since this case 
is studied in the literature also for very large lattice sizes 
Jl2| . Particular attention must be paid in performing 
simulations with dissipation, because if the dissipation is 
too small, £ can become larger than L leading to spu- 
rious results for the cutoff. It is easy to recognize that 
diminishing the dissipation rates is similar to increasing 
the system sizes; in both cases the average avalanche size 
is increasing. 

Our simulations provide two independent estimates of 
the exponent r by extrapolating the power law behavior 
for different sizes L and finite dissipation rates. The nu- 
merical determination of an overall power law behavior 
determined with a ten percent accuracy is an easy task. 
On the contrary, to increase the accuracy of an order of 
magnitude requires a very careful data treatment. We 
noticed that the individuation of the straight portion of 
the probability distribution is a very delicate point in the 
accurate evaluation of the exponent r. In particular, even 
innocuous smoothing procedures give rise to impressive 
systematic bias. In fact, the fit of the exponent r suffers 
from strong systematic errors due to the lower and upper 
cutoff of the distribution. For this reason, we perform a 
local slope analysis of the raw data by studying the be- 
havior of the logarithmic derivative of each avalanche dis- 
tribution. In this way, it is possible to identify a plateau 
in which the local slope is almost constant. This plateau 
defines the range of s we can use for a meaningful de- 
termination of the exponent r. Naturally, this range is 
increasing for larger sizes L and smaller dissipation rates 
e. Nevertheless, the measurements of r presents strong 
finite size effects especially in d — 2. In this case the 
exponent t seems to suffer from logarithmic corrections 
with the size L; i.e. t(L) = t — const/ log L. In d > 2, 
the numerical evidences show a much faster convergence 
estimated as t(L) = t — constL~ 2 . In the literature, the 
asymptotic estimates of r are obtained through extrapo- 
lation from the previous functional behavior g,[l^,[l2|,[l3| . 
For greater accuracy we used also a new extrapolation 
procedure devised in Ref. This procedure improves 
the determination of the exponent by using the functional 
form of the corrections for the direct determination of r 
by comparing different size samples. In table I we report 
the asymptotic values of the exponent r for 2 < d < 5. 
The values are ingood agreement with previous estimate 
from ref.s |9|, [l0| , [l2|1 . In addition, it appears from the re- 
sults of table I that also in d > 4 the measured value is 
not definitely converged on the MF result. The values 
extrapolated in presence of finite dissipation rates e have 
a small systematic discrepancy with respect to the values 
obtained in the usual extrapolation procedure. However 
this is to be ascribed to the different boundary conditions 
used in the simulations. It is worth to remark that, as 
already pointed out by other authors |ll[ , the sole anal- 
ysis of r can be misleading, since this exponent is not 
very sensible to the variations of the dimension d, as well 
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as variations of universality class The exponent, in 
fact, suffers a maximum variation of around 10 per cent 
with respect to its MF value. The simple analysis of this 
exponent is therefore not always determinant in the dis- 
crimination of many of the crucial properties of sandpilc 
models. 

In order to provide another independent estimate of 
the exponents r, D and a, we perform a data collapse 
analysis, which turns out to be very powerful in this 
case. Under the finite size scaling assumptions, the distri- 
butions P(s,L) and P(s,e) collapse onto a single curve 
if we rescale properly the variables. Thus, by defining 
P qx = P(s,x)/s~ T and q L = sL~ D (q t = se 1/<T ) we ob- 
tain that all data must collapse onto the universal func- 
tion: 

P q * = G(q x ). (3) 

The exponent r controls the rescaling of the vertical axis, 
while the exponents D and a define the rescaling of the 
horizontal axis. Similar universal function can be ob- 
tained by using as rescaling variables L or e, thus obtain- 
ing P(s, x)L Dt = T{sL-°) andP(s,i)r T /" = W(se 1 / CT ) 
The same analysis can be also performed on the inte- 
grated distribution P(s* > s), that is usually less noisy. 
In this case the power law behavior is governed by the 
exponent r — 1. In order to test carefully the numerical 
data we repeated the data collapse analysis by using all 
the previous data collapse forms as well as a direct fitting 
procedure. We show in figure 1 and 2 the data collapse 
for the conservative and dissipative BTW model in d = 4. 
We obtain very precise collapses, which are very sensible 
to the tuning of the the various exponents. The evalu- 
ation of exponents by a direct fit obtains results which 
are in perfect agreement with the data collapse analysis. 
In table II we report the values of the various exponents 
in 2 < d < 5. From the present analysis, we verify that 
= Da ~ 2.0 independently on the dimension. As ex- 
pected also the exponent governing the divergence of the 
average size assumes constantly the value 7 ~ 1. It is a 
striking evidence that while the exponents D and a vary 
from d — 2 to d — 5 more than 30 percent with a clear 
trend to the MF values, their product fluctuates of just 
a few percent. This definitely shows that the dynamics 
of sandpile maintains MF features also in low dimensions 
as shown in ref. Furthermore, the constant value of 
v^ 1 = Da provides an additional consistency check for 
reliability of our results. 

Looking at table II, we see a strong indication that 
MF behavior has not yet set in d = 4. In fact, con- 
trary to some recent numerical results Jl3| , we obtain 
that D ~ 3.5 and that a" 1 ~ 1.7. These values, ob- 
tained by data collapse, are undoubtedly far from the 
MF ones. They are also fully compatible with the expo- 
nent r as measured independently from the extrapolation 
procedure. In fact, 7, a and r have to satisfy the scaling 
relations (77 = 2 — r which is fully consistent with 
the measured values. For these reasons, we are confident 



in excluding that d — 4 is the upper critical dimension of 
the sandpile model. 

In order to provide a further check to the previous 
results we also analyzed the avalanche duration distribu- 
tions. The results that will appear in a forthcoming paper 
[fiJIl , confirms the scenario presented in this letter. Here, 
we only report the results concerning d = 4, which are 
important being discriminating for the upper critical di- 
mension. By using the data collapse described previously 
we measured the dynamical critical exponents t c ~ L z 
and t c ~ e~ A defining the divergence of the characteris- 
tic time t c with respect to the system size and dissipation 
rate, respectively. In d = 4 we obtain z — 1.8 ± 0.1 and 
A = 0.90 ± 0.05, that also in this case are different from 
the MF values zmf — 2 and Amf = 1- This again sup- 
ports the claim that d c > 4. 

The value of the upper critical dimension is a long 
standing theoretical question in the study of sandpile 
models. Several theoretical estimates (none of them rig- 
orous) give d c = 4 ||, that has been obtained also from 
recent numerical simulations |JL^|. In contrast, other nu- 
merical studies |l6|] and the analogies with dynamical 
percolation led several authors to conjecture d c = 6. 
^From the analysis of our data, that have been obtained 
using the largest lattice sizes ever used, we can definitely 
say that d c > 4. In d = 5 we note discrepancies between 
the values we measure and MF predictions. However, 
because of the relatively small sizes reached in this case, 
we can not exclude that deviations from the MF behav- 
ior are due to finite size effects. In d = 6 we obtain the 
MF values, but the error bars do not permit a reliable 
discussion of the results. 

The main part of the numerical simulations have been 
run on the Kalix parallel computer fl9|| (a Beowulf 
project at Cagliari Physics Department). We thank G. 
Mula for leading the effort toward organizing this com- 
puter facility. The Center for Polymer Studies is sup- 
ported by NSF. 
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1/a 


v = {Da)- 1 
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0.998 ± 0.001 


2.7 ±0.1 


1.30 ±0.05 


0.48 ±0.03 


3 


0.998 ± 0.001 


3.0 ±0.1 


1.50 ±0.05 


0.50 ±0.02 


4 


0.978 ± 0.003 


3.5 ±0.1 


1.72 ±0.05 


0.49 ±0.02 


5 


0.990 ± 0.003 


3.8 ±0.1 


1.88 ±0.05 


0.50 ±0.02 



TABLE II. Values of the critical exponents in different di- 
mensions. The results obtained with data collapse analysis 
and direct fitting procedure are the same and the errors we 
report is the one estimated by the fitting procedure. 
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FIG. 1. Data collapse analysis of the avalanche size dis- 
tribution of the BTW model in d = 4 with e = and for 
lattice sizes L — 48, 80, 112 and 144. The plot is on a double 
logarithmic scale. 
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1.30 ±0.01 
1.25 ±0.03 


1.33 ±0.01 
1.31 ±0.01 


1.45 ± 0.01 
1.43 ±0.01 


1.51 ±0.01 
1.49 ±0.01 



TABLE I. Exponent r and r £ obtained from different 
sizes and dissipation rates extrapolation, respectively. The 
systematic difference is given by the different boundary con- 
ditions used in the numerical simulations. 



-9-0 

CO 

Q_ 



t=1.43 
a=0.58 



-8.0 -3.0 2.0 7.0 12.0 

1/cr 

se 

FIG. 2. Data collapse analysis of the avalanche size dis- 
tribution of the BTW model in d = 4 with L = 144 and for 
dissipation rates e > 10~ 1 .The plot is on a double logarithmic 
scale. 
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